clear; clc; close all;

% 总迭代步数
N = 180;
n = 0:N;
threshold = 1e-8;  % 设定阈值，根据情况调整

%% (a) 绘制 i(n) 和 v(n) 的迭代序列
A = 0.1;          % 幅值
omega = 0.1;      % 角频率
q0 = 0;           % 初始状态

i_seq = A * sin(omega * n);        % 输入电流序列
q_seq = zeros(1, length(n));       % memristor 状态
v_seq = zeros(1, length(n));       % 输出电压
q_seq(1) = q0;

% 迭代更新
for k = 1:length(n) - 1
    v_seq(k) = cos(q_seq(k)) * i_seq(k);
    q_seq(k+1) = q_seq(k) + i_seq(k);
end

subplot(2,2,1);
plot(n, i_seq, 'b', 'LineWidth', 1.5); grid on;
xlabel('n'); ylabel('i(n)');
hold on; 
plot(n, v_seq, 'r', 'LineWidth', 1.5); grid on;
xlabel('n'); ylabel('v(n)');
title('输出电压');


%% (b) 改变频率，绘制频率依赖的滞回回线
A = 0.1; q0 = 0;
omega = [0.1, 0.2, 0.4];

subplot(2,2,2);
hold on;
for w = omega
    i_temp = A * sin(w * n);
    q_temp = zeros(1, N+1);
    v_temp = zeros(1, N+1);
    q_temp(1) = q0;
    for k = 1:N
        v_temp(k) = cos(q_temp(k)) * i_temp(k);
        q_temp(k+1) = q_temp(k) + i_temp(k);
    end

    % 当 i_temp 的绝对值过小（接近零）时，置为 NaN
    idx = abs( v_temp) < threshold;
    i_temp(idx) = NaN;
    v_temp(idx) = NaN;

    plot(i_temp, v_temp, 'LineWidth', 1.5, 'DisplayName',['\omega=',num2str(w)]);
end
xlabel('i(n)'); ylabel('v(n)');
title(' A=0.1, q_0=0 ');
legend('Location','Best'); grid on;


%% (c) 改变幅值，绘制幅值依赖的滞回回线
omega = 0.1; q0 = 0;
amps = [0.05, 0.1, 0.2];

subplot(2,2,3);
hold on;
for A = amps
    i_temp = A * sin(omega * n);
    q_temp = zeros(1, N+1);
    v_temp = zeros(1, N+1);
    q_temp(1) = q0;
    for k = 1:N
        v_temp(k) = cos(q_temp(k)) * i_temp(k);
        q_temp(k+1) = q_temp(k) + i_temp(k);
    end
    % 当 i_temp 的绝对值过小（接近零）时，置为 NaN
    idx = abs( v_temp) < threshold;
    i_temp(idx) = NaN;
    v_temp(idx) = NaN;
    
    plot(i_temp, v_temp, 'LineWidth', 1.5, 'DisplayName', ['A=', num2str(A)]);
end
xlabel('i(n)'); ylabel('v(n)');
title(' \omega=0.1, q0 = 0');
legend('Location','Best'); grid on;


%% (d) 改变初始状态，绘制初始状态依赖的滞回回线
A = 0.1; omega = 0.1;
q0 = [0.5, 0, -0.5];

subplot(2,2,4);
hold on;
for q = q0
    i_temp = A * sin(omega * n);
    q_temp = zeros(1, N+1);
    v_temp = zeros(1, N+1);
    q_temp(1) = q;
    for k = 1:N
        v_temp(k) = cos(q_temp(k)) * i_temp(k);
        q_temp(k+1) = q_temp(k) + i_temp(k);
    end

    % 当 i_temp 的绝对值过小（接近零）时，置为 NaN
    idx = abs( v_temp) < threshold;
    i_temp(idx) = NaN;
    v_temp(idx) = NaN;

    plot(i_temp, v_temp, 'LineWidth', 1.5, ...
         'DisplayName',['q_0=',num2str(q)]);
end
xlabel('i(n)'); ylabel('v(n)');
title(' A=0.1, \omega=0.1 ');
legend('Location','Best'); grid on; hold off;
